rm(list = ls())
require(foreign)
require(stats4)


multinom=function(vec){
return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2])))))
}

#import data
symmetrydata=read.dta("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Expt1.dta")
symmetryframe=data.frame(id=symmetrydata$var1,user_id=symmetrydata$var2,qn=symmetrydata$var3,qid=symmetrydata$var4,order=symmetrydata$var5,bid=symmetrydata$var6,chosen_act=symmetrydata$var9, true_state=symmetrydata$var10)

#reward value
r=10 

#data frame for letters
symframeletters=subset(symmetryframe,user_id<1700)
symframeletters$correct=(symframeletters$true_state<=13 & symframeletters$chosen_act==1)+(symframeletters$true_state>13 & symframeletters$chosen_act==2)
unique(symframeletters$qid)

#data frame for balls
symframeballs=subset(symmetryframe,user_id>1700)
symframeballs$true_state=symframeballs$true_state+6*(symframeballs$qid==3)+4*(symframeballs$qid==4)+2*(symframeballs$qid==5)
symframeballs$correct=(symframeballs$true_state<=10 & symframeballs$chosen_act==1)+(symframeballs$true_state>10 & symframeballs$chosen_act==2)
unique(symframeballs$qid)

#Variables to be stored
acculin=matrix(rep(0,80),nrow=4,ncol=20)
accupow=matrix(rep(0,80),nrow=4,ncol=20)
accugen=matrix(rep(0,80),nrow=4,ncol=20)
accunoneigh=matrix(rep(0,80),nrow=4,ncol=20)
dataaccu=matrix(rep(0,80),nrow=4,ncol=20)

#Qids for each treatment type
lettersqid=c(3,5,7,8)
ballsqid=c(3,4,5,6)

#number of states used by different treatments
Nvec=c(8,12,16,20)

#States for different treatments
states=list()
states[[1]]=c(7:14)
states[[2]]=c(5:16)
states[[3]]=c(3:18)
states[[4]]=c(1:20)


uidlist=unique(symframeballs$user_id)
indioutframe=data.frame(uid=vector('numeric'),likelin=vector('numeric'),likepow=vector('numeric'),likegen=vector('numeric'),likenoneigh=vector('numeric'),likeinat=vector("numeric"))
for(indii in 1:length(uidlist)){
indiframe=subset(symframeballs, user_id==uidlist[indii])

correct=list()
incorrect=list()
# Make correct vecs (needs changing between treatment types)
for(i in 1:4){
correct[[i]]=vector('numeric')
incorrect[[i]]=vector('numeric')
for(j in 1:Nvec[i]){
correct[[i]][j]=sum(indiframe$correct*(indiframe$qid==ballsqid[i])*(indiframe$true_state==states[[i]][j]))
incorrect[[i]][j]=sum((1-indiframe$correct)*(indiframe$qid==ballsqid[i])*(indiframe$true_state==states[[i]][j]))
}
}

multinomlist=list()
multinomvec=vector("numeric")
for(i in 1:4){
multinomlist[[i]]=vector("numeric")
for(j in 1:Nvec[i]){
multinomlist[[i]][j]=multinom(c(correct[[i]][j],incorrect[[i]][j]))
}
multinomvec[i]=sum((multinomlist[[i]]))
}

multinomconst=sum(multinomvec)

######################
#Linear Mutual Information Costs

accuracylin=function(kappaall,kappaloc,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))
return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

likelihoodlin=function(kappaall,kappaloc){
accuvec1=accuracylin(kappaall,kappaloc,1)
accuvec2=accuracylin(kappaall,kappaloc,2)
accuvec3=accuracylin(kappaall,kappaloc,3)
accuvec4=accuracylin(kappaall,kappaloc,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 
if(indii %in% c(10)){
modellin=mle(likelihoodlin,start=list(kappaall=20,kappaloc=20), method="L-BFGS-B",lower=c(0.1,0.1),upper=c(100,1000))
}else{
modellin=mle(likelihoodlin,start=list(kappaall=5.1,kappaloc=5.1), method="L-BFGS-B",lower=c(0.1,0.1),upper=c(100,1000))
}

lllin=logLik(modellin)[1]

#################################################

#power costs
accuracypow=function(kappaall,kappaloc,rho,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))-N*1/N*log(1/N)
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=1/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))-2*0.5*log(0.5)
}
neighborcost[N/2]=(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))-2*0.5*log(0.5)
return(kappaall*overallcost^rho+kappaloc*sum(neighborcost^rho))
}

#Optimization target
negpayoff=function(accu){

return(-1*(r*sum(accu/(N/2))-cost(accu)))

}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2)-c(1:(N/2))/200,negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

likelihoodpow=function(kappaall,kappaloc,rho){
accuvec1=accuracypow(kappaall,kappaloc,rho,1)
accuvec2=accuracypow(kappaall,kappaloc,rho,2)
accuvec3=accuracypow(kappaall,kappaloc,rho,3)
accuvec4=accuracypow(kappaall,kappaloc,rho,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

if(indii %in% c(10)){
modelpow=mle(likelihoodpow,start=list(kappaall=15,kappaloc=15, rho=0.7), method="L-BFGS-B",lower=c(0.1,0.1,0.1),upper=c(100,100,2))
}else{
modelpow=mle(likelihoodpow,start=list(kappaall=5,kappaloc=5, rho=1), method="L-BFGS-B",lower=c(0.1,0.1,0.1),upper=c(100,100,2))
}

llpow=logLik(modelpow)[1]

#################################################
#Gen entropy approach
accuracygen=function(kappaall,kappaloc,q,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
neighborcost=vector('numeric')

if(q==1){

overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))

}else{

if(q==2){
overallcost=kappaall*(-1/N*sum(log(accu/(N/2))+log((1-accu)/(N/2)))-2*log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(-1/2*((accu[j]+accu[j+1])/2*(log(accu[j]/(accu[j]+accu[j+1]))+log(accu[j+1]/(accu[j]+accu[j+1])))+(2-accu[j]-accu[j+1])/2*(log((1-accu[j])/(2-accu[j]-accu[j+1]))+log((1-accu[j+1])/(2-accu[j]-accu[j+1]) ) )) -2*log(2))
}
neighborcost[N/2]=kappaloc*(-1/2*(log(accu[N/2])+log(1-accu[N/2]))-2*log(2))


}else{
overallcost=kappaall*(N^(1-q)/((q-2)*(q-1))*(sum((accu/(N/2))^(2-q)+((1-accu)/(N/2))^(2-q)))-1/((q-2)*(q-1))-log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[j]+accu[j+1])/2*((accu[j]/(accu[j]+accu[j+1]))^(2-q)+(accu[j+1]/(accu[j]+accu[j+1]))^(2-q))+(2-accu[j]-accu[j+1])/2*(((1-accu[j])/(2-accu[j]-accu[j+1]))^(2-q)+((1-accu[j+1])/(2-accu[j]-accu[j+1]))^(2-q)))-1/((q-2)*(q-1))-log(2))

}
neighborcost[N/2]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[N/2])^(2-q)+(1-accu[N/2])^(2-q))-1/((q-2)*(q-1))-log(2))

}

}

return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

likelihoodgen=function(kappaall,kappaloc,q){
accuvec1=accuracygen(kappaall,kappaloc,q,1)
accuvec2=accuracygen(kappaall,kappaloc,q,2)
accuvec3=accuracygen(kappaall,kappaloc,q,3)
accuvec4=accuracygen(kappaall,kappaloc,q,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

modelgen=mle(likelihoodgen,start=list(kappaall=7,kappaloc=6, q=0.8), method="L-BFGS-B",lower=c(0.01,0.01,0.01),upper=c(100,100,10))


llgen=logLik(modelgen)[1]

######################
#Linear Mutual Information Costs no neighborhoods

accuracynoneigh=function(kappaall,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Information cost function
cost=function(accu){
overallcost=kappaall*(accu*log(accu/(N/2))+(1-accu)*log((1-accu)/(N/2)))
return(overallcost)
}

#Optimization target
negpayoff=function(accuscalar){
return(-1*(r*sum(rep(accuscalar,N/2)/(N/2))-cost(accuscalar)))
}

#Find optimal accuracies
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec=rep(optimout$minimum,N/2)
return(accuvec)
}

likelihoodnoneigh=function(kappaall){
accuvec1=accuracynoneigh(kappaall,1)
accuvec2=accuracynoneigh(kappaall,2)
accuvec3=accuracynoneigh(kappaall,3)
accuvec4=accuracynoneigh(kappaall,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 
if(indii %in% c(8,9,12,16,21)){
modelnoneigh=mle(likelihoodnoneigh,start=list(kappaall=5), method="L-BFGS-B",lower=0.1,upper=100)
}else{
if(indii==17){
modelnoneigh=mle(likelihoodnoneigh,start=list(kappaall=3), method="L-BFGS-B",lower=0.11,upper=100)
}else{
modelnoneigh=mle(likelihoodnoneigh,start=list(kappaall=30), method="L-BFGS-B",lower=0.11,upper=100)
}
}

llnoneigh=logLik(modelnoneigh)[1]

#################################
#Inattentive


accuvec1=rep(0.5,Nvec[1]/2)
accuvec2=rep(0.5,Nvec[2]/2)
accuvec3=rep(0.5,Nvec[3]/2)
accuvec4=rep(0.5,Nvec[4]/2)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
llinat=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))

addframe=data.frame(uid=uidlist[indii],likelin=lllin,likepow=llpow,likegen=llgen,likenoneigh=llnoneigh,likeinat=llinat)

indioutframe=rbind(indioutframe,addframe)

print(indii)

}

indioutframe$AIClin=indioutframe$likelin*-2+4
indioutframe$AICpow=indioutframe$likepow*-2+6
indioutframe$AICgen=indioutframe$likegen*-2+6
indioutframe$AICnoneigh=indioutframe$likenoneigh*-2+2
indioutframe$AICinat=indioutframe$likeinat*-2



indioutframe$minAIC=0
for(i in 1:length(indioutframe$AIClin)){
indioutframe$minAIC[i]=min(c(indioutframe$AIClin[i],indioutframe$AICpow[i],indioutframe$AICgen[i],indioutframe$AICinat[i]),indioutframe$AICnoneigh[i])
}

indioutframe$minAIClin=(indioutframe$AIClin==indioutframe$minAIC)
indioutframe$minAICpow=(indioutframe$AICpow==indioutframe$minAIC)
indioutframe$minAICgen=(indioutframe$AICgen==indioutframe$minAIC)
indioutframe$minAICnoneigh=(indioutframe$AICnoneigh==indioutframe$minAIC)
indioutframe$minAICinat=(indioutframe$AICinat==indioutframe$minAIC)

sum(indioutframe$minAIClin)/length(indioutframe$AIClin)
sum(indioutframe$minAICpow)/length(indioutframe$AIClin)
sum(indioutframe$minAICgen)/length(indioutframe$AIClin)
sum(indioutframe$minAICnoneigh)/length(indioutframe$AIClin)
sum(indioutframe$minAICinat)/length(indioutframe$AIClin)

write.table(indioutframe, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/experiment_4_individuals.csv", sep = ",", col.names = T, append = F)



